1  Prognosemodel Retentie na 1 jaar

ITD | B Communication and Multimedia Design (CMD) - voltijd - versie 1.0

1.1 Methode, data en analyse

1.1.1 Toelichting op de methode

Voor de ontwikkeling van prognosemodellen gebruiken we de aanpak van Tidymodels. Tidymodels is een framework voor het bouwen van een prognosemodel. Hiermee verzekeren we ons van een systematische, herhaalbare en schaalbare aanpak.

1.1.2 Toelichting op de data

De basis voor deze analyse is studiedata van De Haagse Hogeschool (De HHs), verrijkt door het lectoraat LTA. De data bevat informatie over de inschrijvingen van studenten in het eerste jaar van de opleiding:

1. Demografische kenmerken: geslacht, leeftijd, reistijd en SES totaalscore.
2. Vooropleidingskenmerken: toelaatgevende vooropleiding, studiekeuzeprofiel en gemiddeld eindcijfer in de vooropleiding.
3. Aanmeldingskenmerken: aansluiting (direct na diploma, tussenjaar, switch), dag van aanmelding, aantal parallelle studies aan De HHs en collegejaar.

Deze variabelen zijn gekozen omdat we uit eerder onderzoek weten dat ze voorspellende waarde kunnen hebben voor studiesucces (Bakker, 2022) of omdat ze behoren tot sensitieve kenmerken die in fairness analyse gebruikt worden.

# Read the data dictionary
dfData_dictionary <- Get_Data_Dictionary()

# Show the data dictionary
Get_tblData_Dictionary(dfData_dictionary)

Tabel 1.1: Variabelen en mogelijke waarden

Variabele	Toelichting	Mogelijke waarden
Aanmelding	De dag van de aanmelding voor de studie gerekend vanaf 1 september	-366 - 366
Aansluiting	De manier waarop een student instroomt in de opleiding	2e studie, Direct, Na CD, Onbekend, Overig, Switch Extern, Switch Intern, Tussenjaar
APCG	Of de student ten tijde van het behalen van de toelaatgevende vooropleiding in een armoedeprobleemcumulatiegebied woonde	Ja, Nee, Onbekend
Cijfer_CE_Engels	Het gemiddelde cijfer Engels van het centaal examen van de middelbare school	0-10
Cijfer_CE_Natuurkunde	Het gemiddelde cijfer voor Natuurkunde van het centaal examen van de middelbare school	0-10
Cijfer_CE_Nederlands	Het gemiddelde cijfer Nederlands van het centaal examen van de middelbare school	0-10
Cijfer_CE_VO	Het gemiddelde cijfer voor het centaal examen van de middelbare school	0-10
Cijfer_CE_Wiskunde	Het gemiddelde cijfer voor Wiskunde van het centaal examen van de middelbare school	0-10
Cijfer_SE_VO	Het gemiddelde cijfer voor het schoolexamen van de middelbare school	0-10
Collegejaar	Het collegejaar van de inschrijving	2012-2022
Dubbele_studie	Gegeven of de student meer dan 1 studie volgt	Ja, Nee
Geslacht	Het geslacht van de student	M, V
ID	Uniek nummer per student
Leeftijd	De leeftijd van de student op 1 oktober	16-100
Ranking	De positie die de student heeft in de selectie voor de opleiding (alleen bij de B Huidtherapie)	0-350
Reistijd	De reistijd van de student naar de vestiging van de opleiding op een maandag om 9 uur met het openbaar vervoer vanaf de postcode waarop de student woonde ten tijde van het behalen van de toelaatgevende vooropleiding in minuten	0-120
Retentie	Of de student doorstudeert na het eerste studiejaar	Ja, Nee
SES_Arbeid	De sociaaleconomische status score op arbeid van het CBS op basis van de postcode waarop de student woonde ten tijde van het behalen van de toelaatgevende vooropleiding	-1 - 1
SES_Totaal	De sociaaleconomische status score totaal van het CBS op basis van de postcode waarop de student woonde ten tijde van het behalen van de toelaatgevende vooropleiding	-1 - 1
SES_Welvaart	De sociaaleconomische status score op welvaart van het CBS op basis van de postcode waarop de student woonde ten tijde van het behalen van de toelaatgevende vooropleiding	-1 - 1
Studiekeuzeprofiel	Het studiekeuzeprofiel in het voortgezet onderwijs (havo of vwo) of profiel in het mbo	AHO, ALG, CERT, CM, EA, EM, EM&CM, HB, HO, ICT, MedV, NG, NT, NT&NG, Onbekend, TP, TR, TSL, VNL, VS, ZW
Vooropleiding	De toelaatgevende vooropleiding	BD, CD, HAVO, HO, MBO, Overig, VWO

1.1.3 Toelichting op de analyse

We toetsen in deze analyse Retentie na 1 jaar, voortaan Retentie genoemd.

Retentie is gedefinieerd als ingeschreven staan in dezelfde opleiding in een aansluitend collegejaar. Een wisseling van opleidingsvorm binnen de opleiding, bijvoorbeeld van voltijd in jaar 1 naar duaal in jaar 2, geldt ook als retentie. Uitval is het tegenovergestelde van retentie: niet ingeschreven staan in dezelfde opleiding in een aansluitend collegejaar. Een wisseling van opleidingsvorm binnen de opleiding, bijvoorbeeld van voltijd in jaar 1 naar duaal in jaar 2, geldt niet als uitval.

1.2 Voorbereidingen

1.2.1 Laad de data

We laden een subset in van historische data specifiek voor:

Opleiding: ITD | B Communication and Multimedia Design (Synth) (CMD), voltijd, eerstejaars - Retentie na 1 jaar

## Read the data for this study progamme
if(params$use_synthetic_data) {
  dfOpleiding_inschrijvingen_base <- Get_Studyprogram_Enrollments_Synthetic(
    studytrack = params$opleiding,
    studyform = toupper(params$opleidingsvorm_afkorting)
  ) |> 
    mutate(
      INS_Student_UUID_opleiding_vorm = paste(ID, INS_Opleiding, INS_Opleidingsvorm, sep = "_"),
      INS_Opleidingsnaam_huidig = paste(INS_Opleidingsnaam_huidig, "(Synth.)", sep = " ")
    )
} else {
  dfOpleiding_inschrijvingen_base <- get_lta_studyprogram_enrollments_pin(
    board = "HHs/Inschrijvingen",
    faculty = params$faculteit,
    studyprogram = params$opleidingsnaam,
    studytrack = params$opleiding,
    studyform = toupper(params$opleidingsvorm),
    range = "eerstejaars")
}

## Rearrange the levels
Set_Levels(dfOpleiding_inschrijvingen_base)

dfOpleiding_inschrijvingen_base <- dfOpleiding_inschrijvingen_base |>  
  
  ## Create a simple success variable
  Mutate_Retentie(sSucces_model) |>
  
  ## Convert the success variable into a factor
  mutate(SUC_Retentie = as.factor(SUC_Retentie)) |> 

  ## Special possibly based on the propaedeutic diploma
  # Filter_Propedeusediploma(sPropedeusediploma) |>

  ## Make the Dual Study variable a Yes/No variable
  mutate(INS_Dubbele_studie = ifelse(INS_Aantal_inschrijvingen > 1, "Ja", "Nee")) |>  

  ## Remove INS_Aantal_inschrijvingen
  select(-INS_Aantal_inschrijvingen) |> 

  ## Adjust levels for some variables
  Mutate_Levels(
  c(
    "VOP_Studiekeuzeprofiel_LTA_afkorting",
    "INS_Aansluiting_LTA",
    "VOP_Toelaatgevende_vooropleiding_soort"
  ),
    list(lLevels_skp, lLevels_vop, lLevels_vop)
  )
  
## B Huidtherapie: Filter on only students with a grade number (selection)
if(opleiding == "HDT") {
  dfOpleiding_inschrijvingen_base <- dfOpleiding_inschrijvingen_base |> 
    filter(!is.na(RNK_Rangnummer)) 
} 

1.2.2 Selecteer en inspecteer de data

We selecteren eerst de relevante variabelen. We verwijderen daarbij variabelen die maar 1 waarde hebben. We inspecteren de variabelen in een samenvatting in relatie tot retentie en corrigeren daarbij voor multiple testing; de gecorrigeerde significantie-waarden staan vermeld als q-value. Daarnaast inspecteren we de kwaliteit van de data op missende waarden.

lSelect <- c(
    "INS_Student_UUID_opleiding_vorm",
    "CBS_APCG_tf",
    "DEM_Geslacht",
    "DEM_Leeftijd_1_oktober",
    "GIS_Tijd_fiets_OV",
    "INS_Collegejaar",
    "INS_Dagen_tussen_aanmelding_en_1_september",
    "INS_Dubbele_studie",
    "INS_Aansluiting_LTA",
    "SES_Deelscore_arbeid",
    "SES_Deelscore_welvaart",
    "SES_Totaalscore",
    "SUC_Retentie",
    "VOP_Gemiddeld_cijfer_cijferlijst",
    "VOP_Gemiddeld_eindcijfer_VO_van_de_hoogste_vooropleiding_voor_het_HO",
    "VOP_Cijfer_CE1_nederlands",
    "VOP_Cijfer_CE1_engels",
    "VOP_Cijfer_CE_proxy_wiskunde",
    "VOP_Cijfer_CE1_natuurkunde",
    "VOP_Studiekeuzeprofiel_LTA_afkorting",
    "VOP_Toelaatgevende_vooropleiding_soort"
  )

## B Huidtherapie: add the variable RNK_Rangnummer
if(opleiding == "HDT") {
  lSelect <- c(lSelect, "RNK_Rangnummer")
}

## Create a subset
dfOpleiding_inschrijvingen <- dfOpleiding_inschrijvingen_base |>
  
  ## Select the relevant variables
  select_at(lSelect) |>
  
  ## Rename variables for more readable names
  rename(
    ID                    = INS_Student_UUID_opleiding_vorm,
    Geslacht              = DEM_Geslacht,
    Leeftijd              = DEM_Leeftijd_1_oktober,
    Reistijd              = GIS_Tijd_fiets_OV,
    Dubbele_studie        = INS_Dubbele_studie,
    Collegejaar           = INS_Collegejaar,
    Aanmelding            = INS_Dagen_tussen_aanmelding_en_1_september,
    Aansluiting           = INS_Aansluiting_LTA,
    APCG                  = CBS_APCG_tf,
    SES_Arbeid            = SES_Deelscore_arbeid,
    SES_Welvaart          = SES_Deelscore_welvaart,
    SES_Totaal            = SES_Totaalscore,          
    Retentie              = SUC_Retentie,
    Cijfer_SE_VO          = VOP_Gemiddeld_cijfer_cijferlijst,
    Cijfer_CE_VO          = VOP_Gemiddeld_eindcijfer_VO_van_de_hoogste_vooropleiding_voor_het_HO,
    Cijfer_CE_Nederlands  = VOP_Cijfer_CE1_nederlands,
    Cijfer_CE_Engels      = VOP_Cijfer_CE1_engels,
    Cijfer_CE_Wiskunde    = VOP_Cijfer_CE_proxy_wiskunde,
    Cijfer_CE_Natuurkunde = VOP_Cijfer_CE1_natuurkunde,
    Studiekeuzeprofiel    = VOP_Studiekeuzeprofiel_LTA_afkorting,
    Vooropleiding         = VOP_Toelaatgevende_vooropleiding_soort
  ) |> 
  
  ## Adjust CBS_APCG_tf to a factor
  mutate(APCG = case_when(APCG == TRUE ~ "Ja",
                          APCG == FALSE ~ "Nee",
                          .default = "Onbekend")) |>

  ## Indicate where missing numbers are in VO
  Mutate_Cijfers_VO() |>
  
  ## Remove variables, where there is only 1 value
  select(where(~ n_distinct(.) > 1)) |>
  
  arrange(Collegejaar, ID)

## B Huidtherapie: rename the variable RNK_Rangnummer
if(opleiding == "HDT") {
  dfOpleiding_inschrijvingen <- dfOpleiding_inschrijvingen |> 
    rename(Rangnummer = RNK_Rangnummer)
} 

dfOpleiding_inschrijvingen <- dfOpleiding_inschrijvingen |> 
 ltabase::sort_distinct()

## Remove the base dataset
## rm(dfOpleiding_inschrijvingen_base)

Tabel 1.2: Variabelen in relatie tot de uitkomstmaat: params$succes

	Retentie
Variabele	Ja, N=1002 (62%)1	Nee, N=611 (38%)1	p-value2	q-value3	Totaal, N = 16131
Aanmelding	135,80 (62,44)	112,36 (68,57)	<0,001	<0,001***	126,92 (65,80)
Aansluiting			<0,001	0,012***
2e Studie	6 (40%)	9 (60%)			15 (100%)
Direct	455 (60%)	301 (40%)			756 (100%)
Na CD	14 (74%)	5 (26%)			19 (100%)
Onbekend	0 (NA%)	0 (NA%)			0 (NA%)
Overig	0 (NA%)	0 (NA%)			0 (NA%)
Switch extern	254 (57%)	191 (43%)			445 (100%)
Switch intern	148 (73%)	54 (27%)			202 (100%)
Tussenjaar	125 (71%)	51 (29%)			176 (100%)
APCG			<0,001	0,013***
Ja	224 (63%)	133 (37%)			357 (100%)
Nee	721 (64%)	411 (36%)			1.132 (100%)
Onbekend	57 (46%)	67 (54%)			124 (100%)
Cijfer_CE_Engels	6,93 (1,18)	7,07 (1,10)	0,16	>0,99	6,98 (1,15)
Cijfer_CE_Engels_missing			0,70	>0,99
Ja	533 (62%)	331 (38%)			864 (100%)
Nee	469 (63%)	280 (37%)			749 (100%)
Cijfer_CE_Natuurkunde	6,18 (0,96)	6,31 (0,95)	0,21	>0,99	6,23 (0,95)
Cijfer_CE_Natuurkunde_missing			0,68	>0,99
Ja	867 (62%)	533 (38%)			1.400 (100%)
Nee	135 (63%)	78 (37%)			213 (100%)
Cijfer_CE_Nederlands	5,94 (0,85)	5,96 (0,91)	0,86	>0,99	5,95 (0,87)
Cijfer_CE_Nederlands_missing			0,61	>0,99
Ja	533 (62%)	333 (38%)			866 (100%)
Nee	469 (63%)	278 (37%)			747 (100%)
Cijfer_CE_VO	6,46 (0,37)	6,35 (0,34)	<0,001	<0,001***	6,42 (0,36)
Cijfer_CE_VO_missing			0,59	>0,99
Ja	488 (61%)	306 (39%)			794 (100%)
Nee	514 (63%)	305 (37%)			819 (100%)
Cijfer_CE_Wiskunde	6,35 (1,12)	6,28 (1,07)	0,65	>0,99	6,33 (1,10)
Cijfer_CE_Wiskunde_missing			0,64	>0,99
Ja	544 (62%)	339 (38%)			883 (100%)
Nee	458 (63%)	272 (37%)			730 (100%)
Cijfer_SE_VO	6,46 (0,38)	6,35 (0,36)	<0,001	<0,001***	6,42 (0,38)
Cijfer_SE_VO_missing			0,59	>0,99
Ja	501 (61%)	314 (39%)			815 (100%)
Nee	501 (63%)	297 (37%)			798 (100%)
Dubbele_studie			<0,001	<0,001***
Ja	12 (26%)	34 (74%)			46 (100%)
Nee	990 (63%)	577 (37%)			1.567 (100%)
Geslacht			<0,001	<0,001***
M	521 (57%)	394 (43%)			915 (100%)
V	481 (69%)	217 (31%)			698 (100%)
Leeftijd	20,01 (2,26)	20,19 (2,61)	0,48	>0,99	20,08 (2,40)
Reistijd	37,01 (18,38)	38,84 (22,79)	0,79	>0,99	37,70 (20,17)
SES_Arbeid	0,01 (0,08)	0,01 (0,08)	0,57	>0,99	0,01 (0,08)
SES_Totaal	0,02 (0,28)	0,01 (0,27)	0,45	>0,99	0,02 (0,28)
SES_Welvaart	0,01 (0,13)	0,01 (0,13)	0,36	>0,99	0,01 (0,13)
Studiekeuzeprofiel			0,054	>0,99
AHO	2 (67%)	1 (33%)			3 (100%)
ALG	0 (0%)	1 (100%)			1 (100%)
CERT	1 (100%)	0 (0%)			1 (100%)
CM	78 (72%)	30 (28%)			108 (100%)
EA	36 (67%)	18 (33%)			54 (100%)
EM	118 (56%)	93 (44%)			211 (100%)
EM&CM	152 (70%)	66 (30%)			218 (100%)
HB	8 (62%)	5 (38%)			13 (100%)
HO	27 (57%)	20 (43%)			47 (100%)
ICT	82 (64%)	47 (36%)			129 (100%)
MedV	132 (65%)	72 (35%)			204 (100%)
NG	115 (63%)	68 (37%)			183 (100%)
NT	56 (60%)	38 (40%)			94 (100%)
NT&NG	62 (60%)	41 (40%)			103 (100%)
Onbekend	90 (52%)	83 (48%)			173 (100%)
TP	2 (40%)	3 (60%)			5 (100%)
TR	5 (56%)	4 (44%)			9 (100%)
TSL	4 (80%)	1 (20%)			5 (100%)
VNL	8 (80%)	2 (20%)			10 (100%)
VS	1 (33%)	2 (67%)			3 (100%)
ZW	23 (59%)	16 (41%)			39 (100%)
Vooropleiding			0,010	0,25*
BD	42 (46%)	50 (54%)			92 (100%)
CD	21 (70%)	9 (30%)			30 (100%)
HAVO	548 (64%)	312 (36%)			860 (100%)
HO	27 (53%)	24 (47%)			51 (100%)
MBO	330 (63%)	192 (37%)			522 (100%)
Overig	0 (NA%)	0 (NA%)			0 (NA%)
VWO	34 (59%)	24 (41%)			58 (100%)
1Mean (SD); n (%)
2Wilcoxon rank sum test; Fisher's Exact Test for Count Data with simulated p-value   (based on 2000 replicates); Pearson's Chi-squared test
3*p<0.05; **p<0.01; ***p<0.001

## Load dlookr (temporary)
suppressMessages(library(dlookr))

## Show a summary of the data, sorted by missing values
diagnose(dfOpleiding_inschrijvingen) |> 
  mutate(missing_percent = round(missing_percent, 2),
         unique_rate = round(missing_percent, 2)) |>
  arrange(desc(missing_percent)) |>
  knitr::kable(col.names = c("Variabelen",
                           "Type",
                           "# Missende waarden",
                           "% Missende waarden",
                           "# Unieke waarden",
                           "Ratio unieke waarden"))
## Detach dlookr
detach("package:dlookr", unload = TRUE)

Tabel 1.3: Kwaliteit van de data voor bewerkingen (gesorteerd op missende waarden)

Variabelen	Type	# Missende waarden	% Missende waarden	# Unieke waarden	Ratio unieke waarden
Cijfer_CE_Natuurkunde	numeric	1400	86.79	42	86.79
Cijfer_CE_Wiskunde	numeric	883	54.74	56	54.74
Cijfer_CE_Nederlands	numeric	866	53.69	45	53.69
Cijfer_CE_Engels	numeric	864	53.56	57	53.56
Cijfer_SE_VO	numeric	815	50.53	23	50.53
Cijfer_CE_VO	numeric	794	49.23	23	49.23
Studiekeuzeprofiel	factor	173	10.73	21	10.73
SES_Welvaart	numeric	126	7.81	355	7.81
SES_Arbeid	numeric	125	7.75	261	7.75
SES_Totaal	numeric	125	7.75	459	7.75
Reistijd	numeric	20	1.24	362	1.24
Aanmelding	numeric	0	0.00	266	0.00
Aansluiting	factor	0	0.00	6	0.00
APCG	character	0	0.00	3	0.00
Cijfer_CE_Engels_missing	character	0	0.00	2	0.00
Cijfer_CE_Natuurkunde_missing	character	0	0.00	2	0.00
Cijfer_CE_Nederlands_missing	character	0	0.00	2	0.00
Cijfer_CE_VO_missing	character	0	0.00	2	0.00
Cijfer_CE_Wiskunde_missing	character	0	0.00	2	0.00
Cijfer_SE_VO_missing	character	0	0.00	2	0.00
Collegejaar	numeric	0	0.00	11	0.00
Dubbele_studie	character	0	0.00	2	0.00
Geslacht	factor	0	0.00	2	0.00
ID	character	0	0.00	1613	0.00
Leeftijd	numeric	0	0.00	17	0.00
Retentie	factor	0	0.00	2	0.00
Vooropleiding	factor	0	0.00	6	0.00

1.2.3 Bewerk de data

• Uit de eerste diagnose blijkt dat niet alle variabelen goed genoeg zijn voor het bouwen van een prognosemodel: er zijn missende waarden en niet alle veldtypes zijn geschikt.
• Om bias te voorkomen verwijderen we geen rijen met missende waarden, maar vullen die op (imputatie). We bewerken de data zo dat alle missende waarden worden opgevuld: bij numerieke waarden met het gemiddelde en bij categorische variabelen met ‘Onbekend’.
• We passen het type van sommige variabelen aan, zodat ze in het model gebruikt kunnen worden: tekstvelden zetten we om naar factor (een categorische variabele); logische variabelen (Ja/Nee) zetten we om naar een numerieke variabele (1/0).
• De uitkomstvariabele, Retentie, leiden we af van de variabele SUC_Uitval_aantal_jaar_LTA. Als de waarde daar 1 is, is de student na 1 jaar uitgevallen, 2 na 2 jaar, etc. Zolang de waarde daar 0 is, is de student niet uitgevallen.
• Een fictief studentnummer (INS_Student_UUID_opleiding_vorm) gebruiken we, zodat we - als er afwijkende resultaten zijn - de dataset gericht kunnen onderzoeken als dat nodig is.

## Edit the data
dfOpleiding_inschrijvingen <- dfOpleiding_inschrijvingen |> 
  
  ## Imputate all numeric variables with the mean
  mutate(across(where(is.numeric), ~ ifelse(
    is.na(.x),
    mean(.x, na.rm = T),
    .x
  )) ) |>
  
  ## Convert character variables to factor
  mutate(across(where(is.character), as.factor)) |> 
  
  ## Convert logical variables to 0 or 1
  mutate(across(where(is.logical), as.integer)) |>
  
  ## Fill in factors missing values with “Unknown”
  mutate(across(where(is.factor), ~ suppressWarnings(
    fct_explicit_na(.x, na_level = "Onbekend")
  ))) |> 
  
  ## Rearrange the columns so that Retentie is in front
  select(Retentie, everything()) 

## View the data
## glimpse(dfOpleiding_inschrijvingen) 

## Load dlookr (temporary)
suppressMessages(library(dlookr))

## Diagnose the data
diagnose(dfOpleiding_inschrijvingen) |> 
  mutate(missing_percent = round(missing_percent, 2),
         unique_rate = round(unique_rate, 2)) |>
  knitr::kable(col.names = c("Variabelen",
                           "Type",
                           "# Missende waarden",
                           "% Missende waarden",
                           "# Unieke waarden",
                           "Ratio unieke waarden"))
## Detach dlookr
detach("package:dlookr", unload = TRUE)

Tabel 1.4: Kwaliteit van de data na bewerkingen (gesorteerd op missende waarden)

Variabelen	Type	# Missende waarden	% Missende waarden	# Unieke waarden	Ratio unieke waarden
Retentie	factor	0	0	2	0.00
Aanmelding	numeric	0	0	266	0.16
Aansluiting	factor	0	0	6	0.00
APCG	factor	0	0	3	0.00
Cijfer_CE_Engels	numeric	0	0	57	0.04
Cijfer_CE_Engels_missing	factor	0	0	2	0.00
Cijfer_CE_Natuurkunde	numeric	0	0	42	0.03
Cijfer_CE_Natuurkunde_missing	factor	0	0	2	0.00
Cijfer_CE_Nederlands	numeric	0	0	45	0.03
Cijfer_CE_Nederlands_missing	factor	0	0	2	0.00
Cijfer_CE_VO	numeric	0	0	23	0.01
Cijfer_CE_VO_missing	factor	0	0	2	0.00
Cijfer_CE_Wiskunde	numeric	0	0	56	0.03
Cijfer_CE_Wiskunde_missing	factor	0	0	2	0.00
Cijfer_SE_VO	numeric	0	0	23	0.01
Cijfer_SE_VO_missing	factor	0	0	2	0.00
Collegejaar	numeric	0	0	11	0.01
Dubbele_studie	factor	0	0	2	0.00
Geslacht	factor	0	0	2	0.00
ID	factor	0	0	1613	1.00
Leeftijd	numeric	0	0	17	0.01
Reistijd	numeric	0	0	362	0.22
SES_Arbeid	numeric	0	0	261	0.16
SES_Totaal	numeric	0	0	459	0.28
SES_Welvaart	numeric	0	0	355	0.22
Studiekeuzeprofiel	factor	0	0	21	0.01
Vooropleiding	factor	0	0	6	0.00

1.2.4 Bekijk de onderlinge correlaties

Het is verstandig om voorafgaand aan het bouwen van een model te kijken naar de onderlinge correlaties tussen numerieke variabelen. Dit geeft inzicht in de data en kan helpen bij het maken van keuzes voor het model of de duiding van de uitkomsten.

## Create a plot of the intercorrelations in numerical variables
## Remove columns with a standard deviation of 0
dfOpleiding_inschrijvingen |> 
  select(-Collegejaar) |>
  select(where(is.numeric)) |> 
  select_if(~ sd(.) != 0) |>
  cor() |> 
  corrplot::corrplot(
    order = 'hclust', 
    addrect = 4,
    method = "number",  
    tl.cex = 0.8,       
    tl.col = "black",
    diag = FALSE)

Figuur 1.1: Correlatiematrix

1.2.5 Bouw de trainingset, validatieset en testset

• De data is nu geschikt om een prognosemodel mee te bouwen.
• Om het model te bouwen, testen en valideren, splitsen we de data in drie delen van 60%, 20% en 20%. We doen dit op zo’n manier, dat elk deel ongeveer een gelijk aantal studenten bevat dat doorstudeert (dus niet uitvalt).
• We trainen het model op basis van 60% en valideren de modellen tijdens het trainen op de overige 20% (de validatieset).
• De verdeling van de training- en validatieset muteren we 10x (10 folds) om te voorkomen dat het model te veel leert van de trainingset en daardoor slecht presteert op de validatieset.
• Als het model klaar is, testen we het op de 20% studenten uit de testset. De testset blijft dus de gehele tijd ongemoeid, zodat we overfitting - een te goed model op bekende data, maar slechte presetaties (performance) op onbekende data - voorkomen.
• Een willekeurig, maar vaststaand seed-getal voorkomt dat we bij elke run van het model c.q. deze code een net iets andere uitkomst krijgen.

knitr::include_graphics(here::here("R/images", "voorspelmodel-dataset-lta-hhs.png"))

Figuur 1.2: Splitsing van de dataset in trainingset, validatieset en testset

set.seed(0821)

## Split the data into 3 parts: 60%, 20% and 20%
splits      <- initial_validation_split(dfOpleiding_inschrijvingen,
                                        strata = Retentie,
                                        prop = c(0.6, 0.2))

## Create three sets: a training set, a test set and a validation set
dfRetentie_train      <- training(splits)
dfRetentie_test       <- testing(splits)
dfRetentie_validation <- validation_set(splits)

## Create a resample set based on 10 folds (default)
dfRetentie_resamples  <- vfold_cv(dfRetentie_train, strata = Retentie)

Tabel 1.5: Verhouding van de uitkomstvariabele in de training- en testset

Naam	Retentie	Aantal	Proportie
Trainingset	FALSE	366	37.8%
Trainingset	TRUE	601	62.2%
Testset	FALSE	123	38.0%
Testset	TRUE	201	62.0%

1.3 Model I: Logistische Regressie

• Het eerste model is een logistische regressie met penalized likelihood; we gebruiken de glmnet engine voor het bouwen van het model. Penalized likelihood is een techniek die helpt bij het voorkomen van overfitting: het voegt voor elke extra variabele een strafterm toe om eenvoudige modellen te belonen. Glmnet is een veelgebruikt package voor het bouwen van logistische regressiemodellen.
• We gebruiken de Area under the ROC Curve (AUC/ROC) als performance metric. De ROC-curve (Receiver Operating Characteristic) is een grafiek die de prestaties van een classificatiemodel afbeeldt door de verhouding tussen de true positives (sensitiviteit) en de false positives (aspecficiteit = 1-specificiteit) te plotten bij verschillende drempelwaarden. De oppervlakte onder deze curve, bekend als de AUC (Area Under the Curve), kwantificeert het onderscheidingsvermogen van het model; een AUC van 1 duidt op een perfect onderscheidend vermogen, terwijl een AUC van 0,5 wijst op een model zonder onderscheidend vermogen.

1.3.1 Maak het model

Eerst bouwen we het model.

## Build the model: logistic regression
lr_mod <- 
  logistic_reg(penalty = tune(), mixture = 1) |> 
  set_engine("glmnet")

1.3.2 Maak de recipe

Vervolgens zetten we meerdere stappen in een ‘recipe’:

• We definiëren de student-ID als ID variabele. Daarmee krijgt deze variabele de rol van uniek rij-kenmerk.
• We verwijderen vervolgens de oorspronkelijke student-ID en het collegejaar uit de data, omdat deze verder niet gebruikt moeten worden in het model.
• We converteren factoren naar dummy variabelen: voor elke categorie wordt er een nieuwe logische variabele (Ja/Nee) aangemaakt.
• We verwijderen variabelen die geen waarde toevoegen: variabelen met uitsluitend nullen.
• We normaliseren numerieke variabelen om ze met elkaar te kunnen vergelijken door ze te centreren en schalen.
• Sterk gecorreleerde waarden verwijderen we nu niet, omdat we later in de analyse de eventuele samenhang met andere variabelen in een prognosemodel nog willen kunnen visualiseren.

## Build the recipe: logistic regression
lr_recipe <- 
  recipe(Retentie ~ ., data = dfRetentie_train) |>  
  update_role(ID, new_role = "ID") |>           ## Set the student ID as an ID variable
  step_rm(ID, Collegejaar) |>                   ## Remove ID and college year from the model
  step_dummy(all_nominal_predictors()) |>       ## Create dummy variables from categorical variables
  step_zv(all_predictors()) |>                  ## Remove zero values
  step_normalize(all_numeric_predictors())      ## Center and scale numeric variables

## Show the recipe
tidy(lr_recipe) |> 
  knitr::kable(col.names = c("Nummer", 
                             "Operatie", 
                             "Type",
                             "Getraind",
                             "Sla over",
                             "ID"))

Tabel 1.6: Recipesteps voor logistische regressie

Nummer	Operatie	Type	Getraind	Sla over	ID
1	step	rm	FALSE	FALSE	rm_KV8Lp
2	step	dummy	FALSE	FALSE	dummy_bUMsj
3	step	zv	FALSE	FALSE	zv_beaLb
4	step	normalize	FALSE	FALSE	normalize_3WTmw

De variabelen die nu nog overblijven zijn:

Tabel 1.7: Resterende variabelen voor logistische regressie na bewerkingen

Aanmelding	APCG_Nee	Studiekeuzeprofiel_MedV
Cijfer_CE_Engels	APCG_Onbekend	Studiekeuzeprofiel_NG
Cijfer_CE_Natuurkunde	Cijfer_CE_Engels_missing_Nee	Studiekeuzeprofiel_NT
Cijfer_CE_Nederlands	Cijfer_CE_Natuurkunde_missing_Nee	Studiekeuzeprofiel_NT.NG
Cijfer_CE_VO	Cijfer_CE_Nederlands_missing_Nee	Studiekeuzeprofiel_Onbekend
Cijfer_CE_Wiskunde	Cijfer_CE_VO_missing_Nee	Studiekeuzeprofiel_TP
Cijfer_SE_VO	Cijfer_CE_Wiskunde_missing_Nee	Studiekeuzeprofiel_TR
Leeftijd	Cijfer_SE_VO_missing_Nee	Studiekeuzeprofiel_TSL
Reistijd	Dubbele_studie_Nee	Studiekeuzeprofiel_VNL
SES_Arbeid	Geslacht_V	Studiekeuzeprofiel_VS
SES_Totaal	Studiekeuzeprofiel_ALG	Studiekeuzeprofiel_ZW
SES_Welvaart	Studiekeuzeprofiel_CM	Vooropleiding_CD
Retentie	Studiekeuzeprofiel_EA	Vooropleiding_HAVO
Aansluiting_Direct	Studiekeuzeprofiel_EM	Vooropleiding_HO
Aansluiting_Na.CD	Studiekeuzeprofiel_EM.CM	Vooropleiding_MBO
Aansluiting_Switch.extern	Studiekeuzeprofiel_HB	Vooropleiding_VWO
Aansluiting_Switch.intern	Studiekeuzeprofiel_HO
Aansluiting_Tussenjaar	Studiekeuzeprofiel_ICT

1.3.3 Maak de workflow

Voor de uitvoering bouwen we een workflow. Daaraan voegen we het model en de bewerkingen in de recipe toe.

## Create the workflow: logistic regression
lr_workflow <- 
  workflow() |>         ## Create a workflow
  add_model(lr_mod) |>  ## Add the model
  add_recipe(lr_recipe) ## Add the recipe

## Show workflow
lr_workflow

══ Workflow ════════════════════════════════════════════════════════════════════
Preprocessor: Recipe
Model: logistic_reg()

── Preprocessor ────────────────────────────────────────────────────────────────
4 Recipe Steps

• step_rm()
• step_dummy()
• step_zv()
• step_normalize()

── Model ───────────────────────────────────────────────────────────────────────
Logistic Regression Model Specification (classification)

Main Arguments:
  penalty = tune()
  mixture = 1

Computational engine: glmnet 

1.3.4 Tune en train het model

Het model moet getuned worden. Dit houdt in dat we de beste parameters voor het model moeten vinden. We maken een grid met verschillende penalty waarden. Daarmee kunnen we vervolgens het beste model selecteren met de hoogste ROC/AUC. We plotten de resultaten van de tuning, zodat we hieruit het beste model kunnen kiezen.

## Create a grid: logistic regression
lr_reg_grid <- tibble(penalty = 10 ^ seq(-4, -1, length.out = 30))

## Train and tune the model: logistic regression
lr_res <- 
  lr_workflow |> 
  tune_grid(dfRetentie_validation,
            grid = lr_reg_grid,
            control = control_grid(save_pred = TRUE),
            metrics = metric_set(roc_auc))

## Plot the results + a red vertical line for the max AUC
lr_plot <- 
  lr_res |> 
  collect_metrics() |> 
  ggplot(aes(x = penalty, y = mean)) + 
  geom_point() + 
  geom_line() + 
  
  ## Make the scale of the x-axis logarithmic
  scale_x_log10(labels = scales::label_number()) +
    theme(
      axis.title.x = element_text(margin = margin(t = 20))
    ) +
  
  # Define the title, subtitle and caption
  labs(
    caption = sCaption,
    x = "Area under the ROC Curve",
    y = "Penalty"
  )
  
  ## Add LTA elements
  lr_plot <- Add_LTA_Theme_Elements(lr_plot, title_subtitle = FALSE)
  
# Find the penalty value with the max AUC
max_auc_penalty <- lr_res |> 
  collect_metrics() |> 
  filter(mean == max(mean)) |> 
  pull(penalty)

# Add the red vertical line to lr_plot
lr_plot_plus <- lr_plot + 
  geom_vline(xintercept = max_auc_penalty, color = "red")

# Find a mean for the max AUC that is higher
max_auc_mean <- lr_res |> 
  collect_metrics() |> 
  filter(mean == max(mean)) |> 
  pull(penalty)

## Print the final plot
lr_plot_plus

Figuur 1.3: Tuning resultaten logistische regressie

1.3.5 Kies het beste model

De prestaties van een model gevisualiseerd met behulp van een ROC curve. De sensitiviteit (True Positive Rate) en specificiteit (True Negative Rate) worden hierin tegenover elkaar uitgezet. De Area under the ROC Curve (AUC/ROC) geeft de prestaties van het model weer. Het model scoort beter naarmate de AUC/ROC dichter bij de 1 ligt, de linker bovenhoek. De linker bovenhoek houdt in dat alle prognoses exact overeenstemmen met de werkelijkheid. Een AUC/ROC van 0,5 betekent dat het model niet beter presteert dan een willekeurige voorspelling.

We gebruiken modellen met een zo hoog mogelijke Area under the ROC Curve (AUC/ROC) en een zo laag mogelijke penalty. Zo kunnen we uit de resultaten het beste model kiezen en visualiseren.

## Show the best model
top_models <-
  lr_res |> 
  show_best(metric = "roc_auc", n = 10) |> 
  mutate(mean = round(mean, 6)) |>
  arrange(penalty) 

top_models|> 
  knitr::kable(col.names = c("Penalty", 
                             "Metriek", 
                             "Estimator",
                             "Gemiddelde",
                             "Aantal",
                             "SE",
                             "Configuratie"))

Tabel 1.8: Model performance voor logistische regressie

Penalty	Metriek	Estimator	Gemiddelde	Aantal	SE	Configuratie
0.0017433	roc_auc	binary	0.668279	1	NA	Preprocessor1_Model13
0.0022122	roc_auc	binary	0.670205	1	NA	Preprocessor1_Model14
0.0028072	roc_auc	binary	0.670820	1	NA	Preprocessor1_Model15
0.0035622	roc_auc	binary	0.672172	1	NA	Preprocessor1_Model16
0.0045204	roc_auc	binary	0.673934	1	NA	Preprocessor1_Model17
0.0057362	roc_auc	binary	0.675369	1	NA	Preprocessor1_Model18
0.0072790	roc_auc	binary	0.676967	1	NA	Preprocessor1_Model19
0.0092367	roc_auc	binary	0.678934	1	NA	Preprocessor1_Model20
0.0117210	roc_auc	binary	0.677336	1	NA	Preprocessor1_Model21
0.0148735	roc_auc	binary	0.669262	1	NA	Preprocessor1_Model22

## Select the best model: logistic regression
lr_best <- 
  lr_res |> 
  collect_metrics() |> 
  filter(mean == max(mean)) |>
  slice(1) 

lr_best|> 
  mutate(mean = round(mean, 6)) |>
  knitr::kable(col.names = c("Penalty", 
                             "Metriek", 
                             "Estimator",
                             "Gemiddelde",
                             "Aantal",
                             "SE",
                             "Configuratie"))

Tabel 1.9: Hoogste model performance voor logistische regressie

Penalty	Metriek	Estimator	Gemiddelde	Aantal	SE	Configuratie
0.0092367	roc_auc	binary	0.678934	1	NA	Preprocessor1_Model20

## Collect the predictions and evaluate the model (AUC/ROC): logistic regression
lr_auc <- 
  lr_res |> 
  collect_predictions(parameters = lr_best) |> 
  roc_curve(Retentie, .pred_FALSE) |> 
  mutate(model = "Logistisch Regressie")

## Plot the ROC curve
Get_ROC_Plot(lr_auc, position = 1)

Figuur 1.4: ROC curve voor logistische regressie

## Determine the AUC of the best model
lr_auc_highest   <-
  lr_res |>
  collect_predictions(parameters = lr_best) |> 
  roc_auc(Retentie, .pred_FALSE)

## Add model name and AUC dfModel_results
dfModel_results <- 
  dfModel_results |>
  add_row(model = "Logistic Regression", auc = lr_auc_highest$.estimate)

1.4 Model II: Tree-based ensemble

• Het tweede model is een random forest: een ensemble van beslisbomen (decision trees). Het is een krachtig model dat goed om kan gaan met complexe data en veel variabelen.
• We gebruiken de ranger engine voor het bouwen van het model.

1.4.1 Bepaal het aantal PC-cores

Omdat een random forest model veel berekeningen vereist, willen we daarvoor alle computerkracht gebruiken die beschikbaar is. Het aantal CPU’s (cores), wat verschilt per computer, bepaalt hoe snel het model getraind kan worden. We bepalen het aantal cores en gebruiken dat bij het bouwen van het model.

## Determine the number of cores
cores <- parallel::detectCores()

1.4.2 Maak het model

We bouwen eerst het model. We gebruiken de rand_forest functie om het model te bouwen. We tunen de mtry en min_n parameters. De mtry parameter bepaalt het aantal variabelen dat per boom wordt gebruikt. De min_n parameter bepaalt het minimum aantal observaties dat in een blad van de boom moet zitten. De functie tune() is hier nog een placeholder om de beste waarden voor deze parameters - die we later bepalen - in te kunnen stellen. We gebruiken 1.000 bomen c.q. versies van het model.

## Build the model: random forest

rf_mod <- 
  rand_forest(mtry = tune(), min_n = tune(), trees = 1000) |> 
  set_engine("ranger", num.threads = cores) |> 
  set_mode("classification")

1.4.3 Maak de recipe

We maken een recipe voor het random forest model. We verwijderen de student ID en het collegejaar uit de data, omdat deze niet moet worden gebruikt in het model. Overige stappen zijn bij een random forest minder relevant in tegenstelling tot een regressiemodel.

## Create the recipe: random forest
rf_recipe <- 
  recipe(Retentie ~ ., data = dfRetentie_train) |> 
  step_rm(ID, Collegejaar)                      ## Verwijder ID en Collegejaar uit het model

## Show the recipe
tidy(rf_recipe) |> 
  knitr::kable(col.names = c("Nummer", 
                             "Operatie", 
                             "Type",
                             "Getraind",
                             "Sla over",
                             "ID"))

Tabel 1.10: Recipesteps voor random forest

Nummer	Operatie	Type	Getraind	Sla over	ID
1	step	rm	FALSE	FALSE	rm_HCcZp

1.4.4 Maak de workflow

We voegen het model en de recipe toe aan de workflow voor dit model.

## Create the workflow: random forest
rf_workflow <- 
  workflow() |> 
  add_model(rf_mod) |> 
  add_recipe(rf_recipe)

## Show workflow
rf_workflow

══ Workflow ════════════════════════════════════════════════════════════════════
Preprocessor: Recipe
Model: rand_forest()

── Preprocessor ────────────────────────────────────────────────────────────────
1 Recipe Step

• step_rm()

── Model ───────────────────────────────────────────────────────────────────────
Random Forest Model Specification (classification)

Main Arguments:
  mtry = tune()
  trees = 1000
  min_n = tune()

Engine-Specific Arguments:
  num.threads = cores

Computational engine: ranger 

1.4.5 Tune en train het model

We trainen en tunen het model in de workflow. We maken een grid met verschillende waarden voor de parameters mtry en min_n. We gebruiken de Area under the ROC Curve (AUC/ROC) als performance metric. Met de resultaten van de tuning kiezen we het beste model.

## Show the parameters that can be tuned
rf_mod

Random Forest Model Specification (classification)

Main Arguments:
  mtry = tune()
  trees = 1000
  min_n = tune()

Engine-Specific Arguments:
  num.threads = cores

Computational engine: ranger 

## Extract the parameters being tuned
extract_parameter_set_dials(rf_mod)

Collection of 2 parameters for tuning

 identifier  type    object
       mtry  mtry nparam[?]
      min_n min_n nparam[+]

Model parameters needing finalization:
   # Randomly Selected Predictors ('mtry')

See `?dials::finalize` or `?dials::update.parameters` for more information.

## Determine the seed
set.seed(2904)

## Build the grid: random forest
rf_res <- 
  rf_workflow |> 
  tune_grid(dfRetentie_validation,
            grid = 25,
            control = control_grid(save_pred = TRUE),
            metrics = metric_set(roc_auc))

i Creating pre-processing data to finalize unknown parameter: mtry

1.4.6 Kies het beste model

We evalueren de beste modellen en maken een ROC curve om de performance van het model te visualiseren. Vervolgens vergelijken we de prestaties van de modellen en kiezen daaruit het beste model.

## Show the best models
rf_res |> 
  show_best(metric = "roc_auc", n = 15) |> 
  mutate(mean = round(mean, 6)) |>
  knitr::kable(col.names = c("Mtry", 
                             "Min. aantal", 
                             "Metriek",
                             "Estimator",
                             "Gemiddelde",
                             "Aantal",
                             "SE",
                             "Configuratie"))

Tabel 1.11: Model performance voor random forest

Mtry	Min. aantal	Metriek	Estimator	Gemiddelde	Aantal	SE	Configuratie
2	6	roc_auc	binary	0.700041	1	NA	Preprocessor1_Model22
11	29	roc_auc	binary	0.699836	1	NA	Preprocessor1_Model02
4	26	roc_auc	binary	0.699385	1	NA	Preprocessor1_Model23
5	36	roc_auc	binary	0.699344	1	NA	Preprocessor1_Model09
7	35	roc_auc	binary	0.698074	1	NA	Preprocessor1_Model14
8	31	roc_auc	binary	0.698033	1	NA	Preprocessor1_Model24
11	21	roc_auc	binary	0.697951	1	NA	Preprocessor1_Model15
2	32	roc_auc	binary	0.696434	1	NA	Preprocessor1_Model03
3	8	roc_auc	binary	0.695492	1	NA	Preprocessor1_Model08
6	23	roc_auc	binary	0.695246	1	NA	Preprocessor1_Model17
13	25	roc_auc	binary	0.694262	1	NA	Preprocessor1_Model25
22	38	roc_auc	binary	0.692295	1	NA	Preprocessor1_Model01
16	34	roc_auc	binary	0.691926	1	NA	Preprocessor1_Model07
18	40	roc_auc	binary	0.690369	1	NA	Preprocessor1_Model16
15	20	roc_auc	binary	0.690287	1	NA	Preprocessor1_Model20

## Plot the results
autoplot <- autoplot(rf_res) +
  theme_minimal() +
  labs(
    y = "roc/auc",
    caption = sCaption
  )
  
  ## Add LTA elements
  autoplot <- Add_LTA_Theme_Elements(autoplot, title_subtitle = FALSE)
  
  print(autoplot)

Figuur 1.5: Model performance random forest

## Select the best model
rf_best <- 
  rf_res |> 
  select_best(metric = "roc_auc")

rf_best|> 
  knitr::kable(col.names = c("Mtry", 
                             "Min. aantal", 
                             "Configuratie"))

Tabel 1.12: Hoogste model performance voor random forest

Mtry	Min. aantal	Configuratie
2	6	Preprocessor1_Model22

## Collect the predictions
rf_res |> 
  collect_predictions() |> 
  head(10) |>
  mutate(.pred_FALSE = scales::percent(.pred_FALSE, accuracy = 0.1),
         .pred_TRUE = scales::percent(.pred_TRUE, accuracy = 0.1)) |>
  knitr::kable(col.names = c("% Voorsp. FALSE", 
                             "% Voorsp. TRUE", 
                             "ID",
                             "Rij",
                             "Mtry", 
                             "Min. aantal", 
                             "Retentie",
                             "Configuratie"))

Tabel 1.13: Predicties voor random forest

% Voorsp. FALSE	% Voorsp. TRUE	ID	Rij	Mtry	Min. aantal	Retentie	Configuratie
53.7%	46.3%	validation	968	22	38	TRUE	Preprocessor1_Model01
45.4%	54.6%	validation	969	22	38	TRUE	Preprocessor1_Model01
41.4%	58.6%	validation	970	22	38	TRUE	Preprocessor1_Model01
41.6%	58.4%	validation	971	22	38	TRUE	Preprocessor1_Model01
33.4%	66.6%	validation	972	22	38	TRUE	Preprocessor1_Model01
51.1%	48.9%	validation	973	22	38	TRUE	Preprocessor1_Model01
21.6%	78.4%	validation	974	22	38	FALSE	Preprocessor1_Model01
50.0%	50.0%	validation	975	22	38	TRUE	Preprocessor1_Model01
48.1%	51.9%	validation	976	22	38	FALSE	Preprocessor1_Model01
22.4%	77.6%	validation	977	22	38	TRUE	Preprocessor1_Model01

## Determine the AUC/ROC curve
rf_auc <- 
  rf_res |> 
  collect_predictions(parameters = rf_best) |> 
  roc_curve(Retentie, .pred_FALSE) |> 
  mutate(model = "Random Forest")

## Plot the ROC curve
Get_ROC_Plot(rf_auc, position = 2)

## Determine the AUC of the best model
rf_auc_highest   <-
  rf_res |>
  collect_predictions(parameters = rf_best) |> 
  roc_auc(Retentie, .pred_FALSE)

## Add model name and AUC to dfModel_results
dfModel_results <- 
  dfModel_results |>
  add_row(model = "Random Forest", 
          auc = rf_auc_highest$.estimate)

Figuur 1.6: ROC curve voor random forest

1.5 De uiteindelijke fit

• In de laatste stap van deze analyse maken we het model definitief.
• We testen het model op de testset en evalueren het model met metrieken en de Variable Importance (VI). De VI kwantificeert de bijdrage van elke variabele aan de voorspellende kracht van een model. Het identificeert welke variabelen significant zijn voor de modelprestaties, wat essentieel is voor het interpreteren en optimaliseren van een model (Van der Laan, 2006). Methoden zoals de Shapley-waarde en permutation importance worden vaak toegepast om dit belang te meten. Op deze methoden komen we terug in het volgende hoofdstuk.

1.5.1 Combineer de AUC/ROC curves en kies het beste model

Eerst combineren we de AUC/ROC curves van de modellen om ze te vergelijken. We kiezen het beste model op basis van de hoogste AUC/ROC.

## Combine the AUC/ROC curves to compare the models
Get_ROC_Plot(list(lr_auc, rf_auc))

Figuur 1.7: Gecombineerde ROC curves

## Determine which of the models is best based on highest AUC/ROC
dfModel_results <- dfModel_results |>
  mutate(number = row_number()) |> 
  mutate(best = ifelse(auc == max(auc), TRUE, FALSE)) |> 
  arrange(number)

## Determine the best model
sBest_model     <- dfModel_results$model[dfModel_results$best == TRUE]
sBest_model_auc <- round(dfModel_results$auc[dfModel_results$best == TRUE], 4)

Het beste model is het Random Forest model met een AUC/ROC van 0.7. We ronden de analyse verder af met het Random Forest model.

1.5.2 Maak het finale model

We maken het finale model op basis van de beste parameters die we hebben gevonden. Door in de engine bij importance de impurity op te geven, wordt het beste random forest model gekozen om de data definitief mee te classificeren.

## Test the developed model on the test set
## Determine the optimal parameters

## Build the final models
last_lr_mod <-
    logistic_reg(penalty = lr_best$penalty,
                 mixture = 1) |>
    set_engine("glmnet") |>
    set_mode("classification")

last_rf_mod <-
    rand_forest(mtry = rf_best$mtry,
                min_n = rf_best$min_n,
                trees = 1000) |>
    set_engine("ranger", num.threads = cores, importance = "impurity") |>
    set_mode("classification")

1.5.3 Maak de workflow

We voegen het model toe aan de workflow en updaten de workflow met het finale model.

## Update the workflows
 last_lr_workflow <- 
    lr_workflow |> 
    update_model(last_lr_mod)

 last_rf_workflow <- 
    rf_workflow |> 
    update_model(last_rf_mod)

1.5.4 Fit het finale model

We voeren de finale fit uit. De functie last_fit past het model toe op de validatieset.

## Perform the final fit
set.seed(2904)

## Make a final fit for both models so we can save it for later use
last_fit_lr <- 
    last_lr_workflow |> 
    last_fit(splits)

last_fit_rf <- 
    last_rf_workflow |> 
    last_fit(splits)

lLast_fits <- list(last_fit_lr, last_fit_rf) |> 
  set_names(c("Logistic Regression", "Random Forest"))

## Determine which model is best
if(sBest_model == "Logistic Regression") {
  last_fit <- last_fit_lr
} else if(sBest_model == "Random Forest") {
  last_fit <- last_fit_rf
}

## Keep results, model results and associated data
sFittedmodels_outputpath <- Get_Model_Outputpath(mode = "last-fits")
saveRDS(lLast_fits, file = sFittedmodels_outputpath)

sModelresults_outputpath <- Get_Model_Outputpath(mode = "modelresults")
saveRDS(dfModel_results, file = sModelresults_outputpath)

sData_outputpath <- Get_Model_Outputpath(mode = "data")
saveRDS(dfOpleiding_inschrijvingen, file = sData_outputpath)

1.5.5 Evalueer het finale model: metrieken en variable importance

We evalueren het finale model op basis van 4 metrieken: 1) accuraatheid, 2) ROC/AUC en 3) de Brier score (de Mean Squared Error) en 4) de Variable Importance (VI). Uit de VI is op te maken welke variabelen het meest bijdragen aan de voorspelling van de uitkomstvariabele.

## Collect the metrics
last_fit |> 
  collect_metrics() |> 
  mutate(.estimate = round(.estimate, 4)) |>
  knitr::kable(col.names = c("Metriek", 
                             "Estimator",
                             "Estimate",
                             "Configuratie"))

Metriek	Estimator	Estimate	Configuratie
accuracy	binary	0.6420	Preprocessor1_Model1
roc_auc	binary	0.6491	Preprocessor1_Model1
brier_class	binary	0.2219	Preprocessor1_Model1

# Extract the variable importance
dfVi <- last_fit |>
  extract_fit_parsnip() |>
  vip::vi() |> 
  arrange(desc(Importance)) |>
  head(20)
  
# Create the plot with fill on the variable 'Importance'
importance_plot <- dfVi |> 
  ggplot(aes(x = reorder(Variable, Importance), 
             y = Importance, 
             fill = Importance)) +
  geom_col(show.legend = FALSE) +
  
  ## Create the title and caption
  labs(x = NULL,
       y = "VI-score",
       caption = sCaption) +
  
  theme_minimal() +
  Set_LTA_Theme() +
  
  theme(
    axis.title.x = element_text(margin = margin(t = 20))
  ) +
  
  coord_flip()
  
  ## Add LTA elements
  importance_plot <- Add_LTA_Theme_Elements(importance_plot, title_subtitle = TRUE)

# Show the plot
print(importance_plot)

Figuur 1.8: Meest voorspellende factoren op basis van de Variable Importance (VI)

1.5.6 Plot de ROC curve

Tot slot visualiseren we de prestaties weer met een ROC curve van het beste model.

## Show the roc curve
auc_lf <- last_fit |> 
  collect_predictions() |> 
  roc_curve(Retentie, .pred_FALSE) |> 
  mutate(model = "Last fit")

Get_ROC_Plot(auc_lf, position = 3)

Figuur 1.9: ROC curve finale model

1.6 Conclusies

1.6.1 Het beste prognosemodel voor deze opleiding

Het beste prognosemodel blijkt het Random Forest model te zijn.

• Van de prognosemodellen die we hebben ontwikkeld om retentie na 1 jaar te voorspellen, had het Random Forest model de hoogste AUC/ROC waarde (0.7).

1.6.2 Mate van accuraatheid en lift

Een prognosemodel moet minimaal beter presteren dan een base-model om waarde op basis van accuraatheid toe te voegen. Het base-model neemt als basis de grootste klasse van de gemiddelde retentie na 1 jaar van de afgelopen jaren. Stel we zouden tegen alle studenten zeggen dat ze hun studie gaan halen, dan is de mate van accuratesse gelijk aan dit base-model. Dit base-model is dus altijd hoger dan de 50% lijn van de AUC/ROC curve, tenzij het base-model toevallig precies 50% is.

knitr::include_graphics(here::here("R/images", "basemodel-lift.png"))

Figuur 1.10: Lift afhankelijk van base-model en accuraatheid

De mate van accuraatheid van het prognosemodel is vrij laag (64.2%).

• Base-model: 62.12% – Voor deze opleiding berekenen we het base-model als volgt. Van alle studenten studeerde 62.12% door; 100% - 62.12% = 37.88% studeerde niet door. De grootste klasse van deze twee, 62.12%, is daarmee de accuratesse van het base-model.
• Accuratesse prognose: 64.2% – Het model voorspelt Retentie na 1 jaar met een accuratesse van 64.2%.
• Lift: 2.08% – Het model scoort in de huidige opbouw met een verschil van 2.08% (de lift) iets beter dan de accuraatheid van het base-model.

1.6.3 Confusion Matrix

## Determine the confusion matrix
confusion_matrix <- last_fit |>
  collect_predictions() |>
  conf_mat(truth = Retentie, estimate = .pred_class) 

dfConf_matrix <- as_tibble(confusion_matrix$table) |>
  rename(Werkelijkheid = Truth) |>
  mutate(Werkelijkheid = ifelse(Werkelijkheid == "TRUE", "Retentie", "Geen retentie"),
         Prediction    = ifelse(Prediction == "TRUE", "Retentie", "Geen retentie"))

pTP  <- Change_Number_Marks((dfConf_matrix$n[4]/sum(dfConf_matrix$n)*100),1)
pFP  <- Change_Number_Marks((dfConf_matrix$n[2]/sum(dfConf_matrix$n)*100),1)
pTN  <- Change_Number_Marks((dfConf_matrix$n[1]/sum(dfConf_matrix$n)*100),1)
pFN  <- Change_Number_Marks((dfConf_matrix$n[3]/sum(dfConf_matrix$n)*100),1)
pACC <- Change_Number_Marks(Last_fit_Accuracy,1)

De prestaties van het model kunnen we verder uitdrukken in een confusion matrix. Hierin zien we de voorspellingen van het model en de werkelijke uitkomsten. De matrix geeft inzicht in de mate van correcte en incorrecte voorspellingen. Ter illustratie werken we de matrix uit voor een voorspelling waarop een bindend studieadvies (BSA) gebaseerd zou kunnen zijn.

knitr::include_graphics(here::here("R/images", "confusion-matrix-retention-lta-hhs.png"))

Figuur 1.11: Confusion matrix in relatie tot BSA

We passen de confusion matrix nu toe op het model dat als beste naar voren kwam. De accuraatheid van dit model is 64,2%. De accuraatheid van het model berekenen we door de som van de diagonaal te berekenen: het aandeel goed voorspelde uitkomsten, Retentie = Retentie (True Positive) en Geen retentie = Geen retentie (True Negative), af te zetten tegen het totaal aantal voorspellingen: 57,4% + 6,8% = 64,2%. (NB. De weergave in deze confusion matrix is diagonaal gespiegeld vergeleken met het voorbeeld.)

confusion_plot <- plot_confusion_matrix(
    dfConf_matrix,
    target_col = "Werkelijkheid",
    prediction_col = "Prediction",
    counts_col = "n",
    palette = "Blues",
    add_sums = TRUE,
    theme_fn = ggplot2::theme_light,
    sums_settings = sum_tile_settings(
      palette = "Greens",
      label = "Totaal",
      tc_tile_border_color = "black"
    )) +
    
    ## Customize the labels
    labs(
      x = "Werkelijke uitkost",
      y = "Voorspelde uitkomst",
      caption = sCaption
    ) +
    
    Set_LTA_Theme()
  
  ## Add LTA elements
  confusion_plot <- Add_LTA_Theme_Elements(confusion_plot, 
                                           title_subtitle = TRUE)
  
  print(confusion_plot)

Figuur 1.12: Confusion matrix ten opzichte van Retentie na 1 jaar

1.6.4 Uitleggen of verklaren?

Naast de accuraatheid van het model is het ook belangrijk om te weten welke factoren het meest bijdragen aan de voorspelling van retentie na 1 jaar. Daarin gaat de vergelijking met de prestaties van het basemodel mank. Dat model geeft op geen enkele manier aan waarom een student een kans op succes heeft, anders dan - ‘dit is gebruikelijk in deze opleiding’.

Ongeacht de mate van accuraatheid, is het voor onderzoek naar kansengelijkheid essentieel om te weten welke factoren het meest bijdragen aan de voorspelling van retentie na 1 jaar. Het gaat erom dat we het belang van de factoren in de voorspellingen kunnen begrijpen en duiden. Machine Learning is hiervoor uitstekend geschikt, omdat het de mogelijkheid biedt om de belangrijkste factoren en hun invloed te leren kennen (Shmueli, 2010; Shmueli & Koppius, 2011).

1.7 Vervolgstappen: Factoranalyse

De volgende stap (stap 2) is een verdiepende analyse van de mate waarin de factoren die we gevonden hebben van invloed zijn op Retentie na 1 jaar. We kijken naar de rangorde, of ze retentie na 1 jaar verhogen of juist verlagen en hoe stabiel de factoren zijn als we in andere volgordes aan het model toevoegen. Om het concreet te maken zullen we het model toepassen op een aantal fictieve studenten, die we opbouwen uit de meeste voorkomende waarden in deze opleiding. Dit is het onderwerp van analyse 2: de Factoranalyse.

 

Verantwoording

Deze analyse maakt deel uit van het onderzoek naar kansengelijkheid van het lectoraat Learning Technology & Analytics van De Haagse Hogeschool: No Fairness without Awareness | Het rapport is door het lectoraat ontwikkeld in Quarto 1.6.39. | Template versie: 0.9.1.9000

 

Copyright

Dr. Theo Bakker, Lectoraat Learning Technology & Analytics, De Haagse Hogeschool © 2023-2025. Alle rechten voorbehouden.

Bakker, T. (2022). Study Progression and Success of Autistic Students in Higher Education. A Longitudinal, Propensity Score-Weighted Population Study. https://doi.org/10.5463/thesis.1

Shmueli, G. (2010). To Explain or to Predict? Statistical Science, 25(3), 289–310. https://doi.org/10.1214/10-sts330

Shmueli, G., & Koppius, O. (2011). Predictive Analytics in Information Systems Research. MIS Quarterly, 35(3), 553. https://doi.org/10.2307/23042796

Van der Laan, M. J. (2006). Statistical inference for variable importance. The International Journal of Biostatistics, 2(1).